##  data analysis '53 paper 
##  spatial probit models


rm(list = ls())

library(readxl)
library(multiwayvcov)
library(lmtest)
library(Hmisc)
library(spatialprobit)

setwd("~/Dropbox/Current projects/1953/Replication archive")


##  choose signal strength measure to employ
signal <- "lrm_4"

##  should municipalities be dropped/recoded to account for possible jamming?
##  jamming = km10, km12, km15; adjust = "drop", "min"
jamming <- ""
adjust <- ""

load(paste("1953data", signal, adjust, jamming, ".RData", sep = ""))

#load("distance matrix.RData")	# d matrix: distance in km between municipalities



##  rescale predictors to avoid numerical problems
dat$distance_to_berlin <- dat$distance_to_berlin/100
dat$signal <- dat$signal/100


##  run probit model to generate model matrix
probit.out <- glm(protest ~	
							county_capital +
							KVP +
							LandForstwirtschaft +
							Bergbau +
							Metallverarbeitung +
							Chemie +
							HolzKunstmassen +
							Verbrauchsgueter +
							Bauwirtschaft +
							Verkehrswesen +
							Handel +
							Dienstleistung +
							log(population) +
							prop_male +
							density +
							pop_change +
							prop_Land_same39 +
							prop_educ_8 +
							prop_educ_10 +
							prop_educ_12 +
							prop_educ_12plus +
							prop_educ_Fachschule +
							prop_educ_Hochschule +
							prop_evangelisch +
							prop_katholisch +
							prop_glaubenslos +
							living_space +
							turnout32 +
							invalid32 +
							NSDAP32 +
							SPD32 +
							KPD32 +
							Z32 +
							DNVP32 +
							DVP32 +

							signal + I(signal^2) + I(signal^3) +
							distance_to_berlin +
							I(distance_to_berlin^2) +
							I(distance_to_berlin^3) +
							as.factor(size) +
							as.factor(district),
							x = TRUE, y = TRUE, model = TRUE,
							data = dat, family = binomial(link = "probit"))

summary(probit.out)

X <- model.matrix(probit.out)[,-1]
range(apply(X,2,sd))



coor <- as.matrix(dat[,c(6,5)])

##  weights matrices using k = 5, 10, 20 neighbors
W5 <- knn2nb(knearneigh(coor, longlat = TRUE, k = 5))
W5 <- nb2mat(W5, style = "W")
W5 <- Matrix(W5, sparse = TRUE, forceCheck = TRUE)

W10 <- knn2nb(knearneigh(coor, longlat = TRUE, k = 10))
W10 <- nb2mat(W10, style = "W")
W10 <- Matrix(W10, sparse = TRUE, forceCheck = TRUE)

W20 <- knn2nb(knearneigh(coor, longlat = TRUE, k = 20))
W20 <- nb2mat(W20, style = "W")
W20 <- Matrix(W20, sparse = TRUE, forceCheck = TRUE)


##  distance-based weights matrices for d = 20
W.d20 <- dnearneigh(coor, d1 = 0, d2 = 20, longlat = TRUE)
W.d20 <- nb2listw(W.d20, style = "W", zero.policy = TRUE)
W.d20 <- as(W.d20, "CsparseMatrix")
table(apply(W.d20, 1, sum))



##  model with k = 5 nearest neighbors
out.5 <- sarprobit(dat$protest ~ X, W = W5, showProgress = TRUE, ndraw = 20000,
		burn.in = 10000, m = 100, thinning = 20)

summary(out.5)
# plot(out.5)
save(out.5, file = "out5.RData")



##  model with k = 10 nearest neighbors
out.10 <- sarprobit(dat$protest ~ X, W = W10, showProgress = TRUE, ndraw = 20000,
		burn.in = 10000, m = 100, thinning = 20)

summary(out.10)
# plot(out.10)
save(out.10, file = "out10.RData")



##  model with k = 20 nearest neighbors
out.20 <- sarprobit(dat$protest ~ X, W = W20, showProgress = TRUE, ndraw = 20000,
		burn.in = 10000, m = 100, thinning = 20)

summary(out.20)
# plot(out.20)
save(out.20, file = "out20.RData")



##  model with d = 20 km euclidean distance
out.d20 <- sarprobit(dat$protest ~ X, W = W.d20, showProgress = TRUE, ndraw = 20000,
		burn.in = 10000, m = 100, thinning = 20)

summary(out.d20)
# plot(out.d20)
save(out.d20, file = "outd20.RData")



########################
##  create results table
########################
source("hl.R")

table1 <- matrix(NA,44,8)
rownames(table1) <- c(	"intercept",
						"county capital",
						"KVP base",
						"$%$ agriculture and forestry",
						"$%$ mining",
						"$%$ metalworking industry",
						"$%$ chemical industry",
						"$%$ woods and plastics",
						"$%$ consumer goods",
						"$%$ construction",
						"$%$ transportation",
						"$%$ commerce",
						"$%$ services",
						"log(population size)",
						"$%$ population male",
						"population density",
						"$%$ population change 1939--50",
						"$%$ population in same {\\it Land}",
						"$%$ 8 years of education",
						"$%$ 10 years of education",
						"$%$ 12 years of education",
						"$%$ more than 12 years of education",
						"$%$ trade school degree",
						"$%$ university degree",
						"$%$ protestant",
						"$%$ catholic",
						"$%$ agnostic",
						"living space per capita (m$^2$)", 
						"$%$ turnout '32",
						"$%$ invalid votes '32",
						"$%$ NSDAP '32",
						"$%$ SPD '32",
						"$%$ KPD '32",
						"$%$ Zentrum '32",
						"$%$ DNVP '32",
						"$%$ DVP '32",
						"RIAS signal strength",
						"RIAS signal strength$^2$",
						"RIAS signal strength$^3$",
						"distance to East Berlin (100 km)",
						"distance to East Berlin$^2$ (100 km)",
						"distance to East Berlin$^3$ (100 km)",
						"rho",
						"district fixed effects")

table1[44,c(1,3,5,7)] <- "Yes"


load("out5.RData")
load("out10.RData")
load("out20.RData")
load("outd20.RData")

table1[c(1:43),1:2] <- summary(out.5)[c(1:42,63),1:2]
table1[c(1:43),3:4] <- summary(out.10)[c(1:42,63),1:2]
table1[c(1:43),5:6] <- summary(out.20)[c(1:42,63),1:2]
table1[c(1:43),7:8] <- summary(out.d20)[c(1:42,63),1:2]

table1.new <- table1

m <- n <- 1
for(m in 1:nrow(table1))	{
	for(n in 1:ncol(table1))	{

table1.new[m,n] <- format(round(as.numeric(table1[m,n]), digits = 3),
						scientific = FALSE, nsmall = 3)

								}
							}

table1.new[44,c(1,3,5,7)] <- "Yes"
table1.new[44,c(2,4,6,8)] <- ""

latex(table1.new, file = "")





